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1 Introduction 



Consider the compressible Navier-Stokes equations for a viscous, heat conducting 
gas. The governing equations can be found in many books. They are in two spatial 
dimensions 

Conservation of Mass 

¥t p + T x {pu) + h {pv) = ^ (la) 

Conservation of x- Momentum 
d d d 

di^ pu) + ^ pu2 + p ^ + ^ puv ^ = ^ 

d r 2 d d d d d 

Conservation of y- Momentum 
d d d 

d r . d d d r 2 d d 

d- X [p{ dy- U + Yx V)] + d-y [ r { V " d~x U)] 
Conservation of Energy 

§ i E + ^[u(E + p)] + ^[v(E + p)]= (Id) 
d d d d 



The dependent variable are the density /?, the velocity field (f = (-u, and the total 
Energy per unit volume E. The total Energy per unit volume is given by 

E = pe + ip(n 2 + v 2 ) (2) 

where e is the internal energy per unit mass. The thermodynamic variables p and 
e are related to the pressure p through the equation of state 

P = P{pi e ) (equation of state) (3) 
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If we assume that the fluid is a perfect gas, the equation of state state is 

p=( 1 -l)pe = (i-l)[E-~p{u 2 + v 2 )] or T = (4) 

where the constants 7 and R are the ratio of specific heats and the gas constant, 
respectively. T is the temperature. For a perfect gas the specific heat at constant 
volume c v and the specific heat at constant pressure c p are related to 7 and R by 

C = * and c r = ^- (5) 
7 — 1 7 — 1 

We assume in the following that the coefficient of viscosity p is constant. We also 
assumed that the coefficient of bulk viscosity is negligible for the fluid, such that 
the second coefficient of viscosity (i is 

P = -\p ( 6 ) 

Furthermore we assume that Fourier's law for heat transfer holds, with a constant 
coefficient of thermal conductivity k. If the gas is polytropic resp. a perfect gas, 
the internal energy e is related to the temperature T by i.e. 

e = c v T (7) 



For a vanishing viscosity p and thermal conductivity k the Navier-Stokes equations 
reduce to the hyperbolic conservation laws for an ideal compressible gas, also 
denoted as Euler equations. In vector form these equations are 



d d d 

—u(x,y,t) + — f(u(x,y,t)) + — g(u(x, y, t)) = 



(8) 



where the conserved variable and flux functions are given by 



u 



M 

pu 
pv 
\EJ 



f(u) 



/ pu 

2 



pv 
puv 



pv + p 
\v(E+p)J 



0) 



pu + p 
puv 
\u(E + p)J 

The velocity field will be denoted by v = (it, v). The density p and pressure p are 
related to the conserved quantities through the equation of state flU) . 
The Euler equations are a system of hyperbolic conservation laws; i.e. for any value 
u o = (po, PoUq, PqVq, Eo) T with positive density po and positive internal energy eo 
the Jacobian matrix 



nDF(u ) = n x Di(u ) +n y Dg(u ) 



(10) 
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is diagonalizable with real eigenvalues for every unit vector n = (n x , n y ). Therefore 
the structure of a plane-wave 



u(x,y,t) = ip(n- (x,y) - st) 



(11) 



propagating at speed s, is independent of the orientation in space; see for example 
[Lev2002]. The Euler equations are rotational symmetric and Galilean invariant. 
A special class of plane- waves is defined through the plane wave Riemann problem; 
an initial value problem for ([5D with piecewise constant initial data, separated by 
a straight line: 



where xq is a given point in (x,y)-plane and Hq a given unit-direction, u; and u r 

are initial states at time t = t n . 

The integral form of the conservation law ([H]) is 



where n(£) = (n x (^),n y (^)) is the outward pointing unit-normal vector on dD at 
a point = ?/(£)) on where £ is the arclength parametrization of 

the boundary dD. D is an arbitrary bounded convex set with a piecewise smooth 
boundary. u R (x(£), i+; n(£)) denotes the one sided limit in time of the solution to 
the plane-wave Riemann problem (|12p at the cell boundary x(£) in the direction 
n(£); i.e. 



If the solution u is smooth, then we can apply the divergence theorem and obtain 
the differential form ([5]) from the integral form (|13p . If a solution of the Euler 
equations contains discontinuities only the integral form is valid. 
In this paper we consider the stability of the plane-wave Riemann problem (|12p for 
the Euler equations and proof that the plane-wave Riemann problem is unstable 
under perturbations. A relation of this instability to the failing of Godunov's and 
Roe's method reported by Quirk [QUI1994], XU [XU1999], Elling [VE2006], Roe 
[ROE2007] is discussed. 

The results of this paper can be extended to general systems of hyperbolic con- 
servation laws, which have at least one genuinely nonlinear characteristic field and 




u/ for no • (x — xo) < 
u r for < no • (x — xq) 



(12) 




(13) 



u K (x(C), t+; n(f )) := lim u K (x (£), t + e; n(£)) on 3D 



(14) 
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two linear degenerate fields with a double eigenvalue. For clarity we restrict the 
description to the Euler equation for a 7-law gas ((SJ) which serves as a model equa- 
tion for more general systems. 

The outline of this paper is as follows. In section 2 we analyze the stability of 
the plane-wave Riemann problem and prove that the solution is unstable under 
perturbations. Viscosity and heat conduction are taken into account in section 3. 
In section 4 the relation to known nonphysical numerical results (carbuncles) in 
Godunov's and Roe's method is discussed and an explanation for carbuncle in- 
stabilities proposed, which is also consistent with Majda's shock stability analysis 
[MA1983]. The last section contains a short summary and some conclusions. 



2 An Instability in the plane-wave 
Riemann Problem 

In this section we analyze the stability of the plane-wave Riemann problem. Since 
the Euler and Navier-Stokes equations are rotational symmetric, we can assume 
without loss of generality that the plane wave moves in the x-direction and that the 
initial states are separated by the line x(y) = 0. Assume the solution of the plane- 
wave Riemann is given by a single discontinuity, then the initial states u; and u r 
of the Riemann problem satisfy the Rankine-Hugoniot jump condition [CF1948]. 

S[u r - Ui] = f(Ur) - /(U/) (15) 

where s is the shock speed. The normal velocity of the discontinuity can be 
computed from (fT5]) . 

Let us assume that the states u/ and u r are connected by a near stationary 1- 
shock wave with normal shock speed s 1 = s 1 (u/,u r ) < 0. Since the shock is near 
stationary, we have — e s < s , for a small positive number e s . Let u r be a small 
perturbation of the right state. The perturbed plane wave u has the initial data 

u{x,y,0) = < (16) 
[u r for U < x 

and consists of a 1-shock, a contact discontinuity and a 3-wave, which is a weak 
rarefaction or weak shock wave. 

In the following we neglect, without loss of generality, the 3-wave such that u mr = 
u r . A smooth function F exists with 

u r = u; +F(ei,e 2 ,e 3 ;ui) (17) 
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where £2 represents a parameter for the strength of the contact discontinuity and 
£3 is a parameter for the strength of a 3- wave; see [SM1983; Chapter 17]. £\ 
represents a parameter for the strength of a 1-shock. We choose £\ such that 

Ur = u J + F(e 1 ,0,0;ui) (18) 

For a perturbation of £2 with £3 = the pressure p and the x-component of the 
velocity are constant behind the 1-shock. We obtain form the Lax shock conditions 
^i+1/2 j ^ ^2(u r ) where A2(u r ) is the second eigenvalue of the Jacobian A(u r ). For 
a near stationary strong shock wave we can assume that < A2(u r ). Since the 
parameter £2 represents a contact discontinuity, the discontinuity speed s 2 is equal 
to the characteristic speed A2(u r ), which is constant across a contact discontinuity; 
i.e. s 2 = A2(u r ) = A2(u r )- Where u r is given by 

u r = ui + F(e 1 ,£ 2 ,0;uj) (19) 

The y-component of the velocity vi and v r enters the solution of the plane- wave Rie- 
mann problem essentially as a parameter. We can first solve the one dimensional 
Riemann problem ignoring the momentum equation for v and then introduce a 
jump in v at the contact discontinuity to obtain the full plane- wave solution. How- 
ever, vi and v r enter the one dimensional Riemann problem through the pressure 
as a function of the total energy, density and velocity. We assume that v r = O(e) 
and vi = 0(e). For a small perturbation O(e) the change in (jl]) is of order £ 2 . 
Therefore we can assume that the wave structure for the perturbed plane-wave 
Riemann problem 

|uf(0— ,y) for x < , . 

u(x,y,0)= _ ' y (20) 
[u^(0+, y) for < x 

persists up to second order in £, where uf(x,y) = (pi, piui, piv e (x,y), Ei) T and 
u e r (x,y) = (p r ,p r u r ,p r v £ (x,y),E r ) T . 



If v £ (x,y) is discontinuous at x = the perturbation introduces a jump in the 
y-component of the velocity at the contact discontinuity. If v r 7^ v\ the per- 
turbed plane-wave Riemann problem contains a tangential or shear instability. In 
[LL1959;§81] it is proved that such a tangential instability in an incompressible 
non viscous flow is absolutely unstable and may lead to a turbulent flow, and 
it is further mentioned that these instabilities also exists in compressible flows. 
However, we assume in the following that the perturbation v £ (x,y) is a smooth 
function. 

We can consider the parameter £1, £2 in (Tl9|) as functions of y. We assume that 
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the variation of e± in y is small enough such that sj +1 ^ 2 j < < u r still holds. This 
variation can be defined independently of the y-component of the velocity and we 
can assume that v is not affected through this perturbation. The flow for this 
perturbed plane-wave Riemann problem is defined through the two-dimensional 
Euler equation ([5]). A change of the parameter £2 in (|19|) affects only the contact 
discontinuity. Since the pressure and the x-component of the velocity are constant 
across a contact discontinuity, we have behind the 1-shock p(x, y, t) = p r = p. 
u(x,y,t) = u r = u r - Thus the Euler equations reduce behind the 1-shock to: 



V 



d d d d 

—p + u r —p + v—p + p—v = (21a) 
Ot ox Oy oy 

d d d . . 

tt-v + u r —v + v—v = 21b) 
ot ox oy 

d d d d 

—E + u r —E + v—E + E—v = (21c) 
ot ox oy oy 

Using ([5]) to rewrite the total energy in ([5]) as a function of the density, pressure and 
velocities and using the other conservative equations in ©, the energy equation 
may be rewritten as a pressure equation 

d d d d d 

w:P + Ur 7r p + v—p + jp{—u + —v) = (22) 
ot Ox oy ox oy 

which for a constant pressure simply reduce to a divergence condition 

Since the x-component of the velocity is constant behind the 1-Shock (|2"Tj) reduces 
to: 

Ft p+Ur oh p+ w=° (24a) 

Wt v + U ^ v = (24b) 

d d d 

—E + u r —E + v—E = (24c) 
ot ox oy 

where the last equation follows from the equation of state ([!]) and from the first 
two. Therefore the perturbed solution satisfies behind the 1-shock the equations 
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with the initial data 

[ u e (x, v) for x < , 
u(x,y,0+) = _ £ , n (25a) 

(x, yj tor U < x 

with 

u e (x,y) = (p (x,y), p (x,y)u r , p (x,y)vQ(x,y), E (x,y)) T (25b) 
and 

u e (x,y) = {p~o{x,y),h{x,y)ur,Po{x 1 y)vl{x,y),E Q (x,y)) T (25c) 
and 

P / v _ JW(7 - 1) + 5Po(aJ.y)(«? + ^6( a:; >y) 2 ) forx<0 
br/(7 - 1) + 2P0(a?. y) 1/) ) for < x 

where Vq(x, y), po(x,y) and po(x,y) are smooth perturbed initial data. 

Proposition 1: a) If the initial tangential velocity Vq does not depend on y then 
the constant pressure solution of the initial value problem (|24|) . (|25|) is given by: 

[poO^ - u r t,y - v e t) for x < u r t 
p[x,y) = < (26a) 
[p~o{x — u r t,y — v £ t) for < u r t 

and 

v £ = vl{x-u r t) (26b) 
and 

E = Vr/{l - 1) + ip(x, y)(u 2 r + (v e f) (26c) 

b) If the initial tangential velocity Vq depends on y a solution does not exist and 
the pressure is not constant. 

Proof: The solution of the plane-wave Riemann problem has behind the 1-shock, 
i.e. for x > s£ +1 ; 2j -t, a constant pressure p r and a constant x-component of the 
velocity u r . Therefore the Euler equations reduces to (I24p and the solution is 
completely defined through the density p, the tangential velocity v, the constant 
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pressure and the constant x-component of the velocity. We have for x < u r t 

d d d 

Wt P + Ur dx- p + V &y P 

= [-u r + u r ]^Po + [-(vP + t ±vF) - u r t^ + - t^)]^Po 

= - [ ¥t v£ + U 4x V£ + Yy V£] % P ° (2?) 
and similar for x > u r t 
d d d 

= [-u r + u r ]^Po + [~(v £ + t^vF) - u r t^ + v%l - vH^)]^po 

= - [ Wi V£ + U ^x V£ + v£ h v£] %^ (28) 
Furthermore 

Q Q Q 

— v £ + u r — v e = [-u r + u r ]— v £ = (29) 
at ox ox 

If follows form (|27p . f|28j) for an initial tangential velocity Vq = Vq(x) that (|26ap 
and ()26bp are solutions of the initial value problem ([25]) , ([25]) . 
We obtain furthermore for the total energy E for x ^ u r t: 

d d d 

at ox oy 

+ pv£[ ¥t v£ + u ^x v£] = 

Therefore ([25]) defines a solution of the Euler equation (jBJ) in smooth parts of the 
flow. 

The jump condition for a plane- wave moving in the x-direction reduces for a con- 
stant pressure and constant x-component of the velocity to 

s(t)[p(x+,y,i) - p(x-,y,t)] = u r {p(x+,y,t) - p(x-,y,t)] 
s(t)u r [p(x+,y,t) - p(x-,y,t)] = (u r ) 2 [p(x+,y,t) - p r (x-,y,t)] 
s(t)[p(x+,y,t)v £ (x+,t) - p(x-,y,t)v e (x-,t)] = 
u r [p(x+,y,t)v e (x+,t) - p(x-,y,t)v £ (x-,t)} 
s(t)[E(x+,y,t) - E(x-,y,t)} = u r [E(x+,y,t) - E(x-,y,t)} 
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where is is the speed of the discontinuity in the x-direction. Therefore we see that 
the jump conditions are satisfied along the curve x = s(t) = u r t, and (|26|) is a 
weak solution of ([5]). 

Since the solution of initial value problem (|24p , (|25p is unique for smooth velocities 
u r and v e , part a) of the proposition follows. 

If the initial data Vq depend on y the unique solution of (|24bp is for every y given 
by (]26bj) . In this case we obtain for the total energy E for x ^ u r t and a constant 
pressure: 

^E + u r ^E + v^E 
at ox ay 

= 2^ + ^m P + Ur d-x P + V dy- p] 

+ pv£[ ^t v£ + U ^x V£ + v£ Ty V£] 

= pv £ v £ ^-v e (30) 
dy 

Thus the energy equation can only be satisfied if v £ does not depend on y. This 
must also hold for t = 0. If v depends on y we obtain from 



Odd d 

w:P + u r —p + v—p = -'jp—v (31) 
at ox oy Oy 

which cannot be satisfied for a constant pressure solution. This proves part b) and 
completes the proof. □ 

In a constant pressure flow the equations for the conservation of energy follows 
from the momentum and the continuity equations. The Euler equations ([5]) there- 
fore simplify to 

¥t p+u iL p+v iry p = Q (32a) 

Odd 

—u + u—u + v—u = (32b) 
at ox oy 

Odd 

—v + u—v + v—v = (32c) 
at ox oy 

If the velocity field is given by u = u r and v = v(x, t) these equations reduce 
to the equations (|24|) for a perturbed flow behind a shock. This system of partial 
differential equations is unstable under perturbations. Let p, u and v be a constant 
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pressure solution of (|32p and p, u, v infinitesimal perturbations. It is assumed that 

p = p + p u=u+u v=v+v 

is also a smooth constant pressure solution of 0. We obtain by neglecting 
quadratic perturbation terms 

d__d__d_ d d _ d d d 

at ax ay ai ax oy ox oy 

d _ d _ d a _ _ a _ _a_ _ a _ _ a _ 
— ti + u—u + f — u + — u + it— — u + u— m + u—u + t> — u = 
ot ox oy ot ox oy ox oy 

d__d__d_ d d d d d 

— V + U — V + V — V + —V + U — V + V — V + U — V + V — V = 

ot ox oy ot ox oy ox oy 
which can be rewritten as a inhomogeneous hyperbolic system 

d _ d _ d d _ d , , 

TT + P + u TrP + v ^~P = ~ U ^~P~ V TTP (33a 
ot ox oy ox oy 

a _a _<9_ _ d _ _ d _ , . 

7— -u -|- u— — u -\- V——U = u— — u v— — u (33b) 
at ax ay ax oy 

d _ d _ d d _ d 

—v + u—v + v—v = —u—v-v—v (33c) 
at ax oy ox Oy 

The solution of the homogeneous system (f3"3"j) is 

p ft =p°(x -ut,y- vt) (34a) 
n ft =u°(x -ut,y- vt) (34b) 
v h =v°(x -ut,y - vt) (34c) 

where /5°, u° and v° are some arbitrary initial divergence free perturbations. Let 

S?(x,y,t) = -u h -^p-v h -^p (35a) 

S*(x,y,t) = -u h -^u-v h -?-u (35b) 

Sy(x,y,t) = -u h ^-v-v h ^-v (35c) 
ax ay 

By formally applying Duhamel's principle we obtain for the inhomogeneous equa- 



lr The results remain valid if pressure perturbations p = p + p are included 
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tion ([53"j) with the source term ([55]) 



p = p h + J S p (x - u(t - r),y - v(t - r),r)dr 

u = u h + S x (x - u(t - r),y - v(t - r),r)dr 
JO 

ax ay 

■y = -t} h + /" S*(a; - u(t - r),y - v(t - t),t)cLt 
Jo 

d_ 

dy' 



v h -t{u h ^-v + v h ^-v) 

OX 1 ' 



(36a) 



(36b) 



(36c) 



where p = p(x,y,t) , u = u(x,y,t) and v = v(x,y,t). Inserting p, u and v into 
(j32[) we can verify that (|36|) defines a solution of the inhomogeneous hyperbolic 
system ([33]) . We obtained, that for a non constant background flow p, u and v, 
perturbation p, u and v grow linear in time. 

The inhomogeneous hyperbolic system ([33]) can be rewritten in vector form as 
d - d - d - 

— w(x, y, t) + A— w(z, y, i) + B— w(x, y, t) + Cw(x, y, t) = (37) 

where the vector of primitive variables is w = (p, u, v) T and the 3x3 matrices 
A(x,y,t), B(x,y,t) and C(x,y,t) are defined by 

o In 














fv 









A = 









B = 





u 




-H 




V° 









V° 










For a smooth background flow ()37|) is a symmetric hyperbolic system, with matrices 
A(x,y, t), B(x,y,t) and C(x,y,t) depending smoothly on x, y and t. For any 
perturbation w with compact support in (x, y) a energy inequality of the form 

|| w(i) ||< exp(M | t |) || w(0) || (38) 

holds, where the constant M depends only on the magnitude of the symmetric 
part of C and of the first derivative of A with respect to x resp. B with respect to 
y; see [LAX2006; Section 4.3]. This result can also be derived directly from (|36|) . 
Now consider the solution (I36p for a case where the initial background density p 
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is not differentiable at x = 0. For a constant x-component of the velocity u = u r 
and a initial tangential velocity vo = Vq(x), which is enforced through a constant 
pressure, the solution is given by (|26p . In this case the definition of the perturbed 
solution (|36p remains valid, if we interpret the partial derivative of p with respect 
to x at the contact discontinuity as a delta function; i.e. 

d _ 

—p(u r t,y,t) = (p r -pi)5(x-u r t) (39) 

where p\ and p r denote the left and right limits at the contact discontinuity. For 
an initial oscillatory velocity perturbation u h , large density perturbations are en- 
forced. In this case the solution grows with t in a manner unlike what we would 
expect from a hyperbolic equation and no energy inequality of the form (|38[) can be 
derived. We denote such a flow in the following as unstable under perturbations. 
Consider for example a vanishing background tangential velocity v, a vanishing per- 
turbed tangential velocity v and a piecewise constant background density p(x, t) 
which is discontinuous at x = ut, then we obtain from (|33a|) for u = u r , 

d d 

^p + u r —p = -u(p r -pi)6(x-u r t) (40) 

and from (|36a|) for the solution 

p = p h - tu h {p r - pi)5{x - u r t) (41) 

For an initial oscillatory velocity perturbation u h large density perturbations are 
enforced. 

For a constant normal velocity u = u r and a piecewise constant background tangen- 
tial velocity v which depends only on x and is discontinuous at x = u r t, equation 
(|33c|) becomes 

d d 

—v + u r —v = —u(v r — vi)5(x — u r t) (42) 

with a weak solution given by (|36cj) ; i.e. 

v = v h - tu h (v r - vi)5(x - Urt) (43) 

For an initial oscillatory velocity perturbation u h , large tangential velocity per- 
turbations are enforced. Linear advection equation of the form (|4U|) . (|42|) with a 
singular source term are discussed by LeVeque in [LEV2002; Section 16.3.1]. 

We showed that the flow behind a perturbed plane stationary shock line is a 
constant pressure flow governed by (|32|) with u = u r and v = v(x,t). Generally, 
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the flow between the nonlinear 1-wave and 3-wave in the Riemann problem is a 
constant pressure flow region. Since a discontinuity in a constant pressure plane 
wave flow region is a contact discontinuity, we obtained: 

If the solution of the plane-wave Riemann problem contains a contact disconti- 
nuity, then the solution is unstable under perturbations. 

Remark: The same perturbation analysis can be applied to an homentrop (con- 
stant entropy), nearly constant density flow, with a not necessary constant pressure 
and divergence. In this case we obtain for the perturbation equations 

Q Q _ Q Q Q 

ot ox oy ox oy 

d _ d _d d _ d _ 1 d 

—u + u—u + v—u = -u—u-v—u-——p (44b) 

ot ox oy ox oy po ox 

d d _ d d d _ 1 d ,. A x 

— v + u— v + v— v = -u— v - v— v - — — p (44c) 
ot ox oy ox oy po oy 

d d d d d 

—p + u—p + v—p + jp(—u+—v) = (44d) 
ot ox oy ox oy 



where po is the constant background density and (u, v) a divergence free back- 
ground velocity field. For a constant background flow po, uo, vo,po, homentrop 
density perturbations are defined through 

p = p A (x - u t,y - v t,t) (45) 

where pa satisfies a wave equation 




where cq denotes the sound speed of the background flow and r refers to the dif- 
ferentiation with respect to the third argument in pA- Since the solutions of the 
wave equation are stable, density perturbations are stable. 

A initial (acoustic) perturbation pA propagates with sound speed relative to the 
background velocity into the flow. If generated at a shock, these acoustic pertur- 
bation propagate with sound speed relative to adjacent (constant) velocities, into 
the flow on both sides of the discontinuity. In contrast to acoustic perturbations, 
density perturbations in a constant pressure region can only exists on one side of 
a shock line. Since in an ideal gas, any density perturbation in a constant pressure 
region is equivalent to a entropy perturbation, the latter disturbances are also de- 
noted as entropy disturbances. 
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In the analysis of the perturbed Riemann problem, the weak 3-wave was neglected. 
Since the change of the entropy across a shock wave is of third order in £3 in (|17|) . 
the weak 3-wave can be regarded as a (discontinuous) homentrop perturbation of 
a constant state. The density of these acoustic perturbation is defined through 
(|4"B]) and do not cause an unstable flow. 



Note: To derive the Rankine-Hugoniot jump condition (|15|) from the integral form 
(|13f) of the conservation law, it is assumed that the integral 

A x /2 f) 

^-Mi } y } t n )di (47) 
-a x /2 at 

approaches zero for limA x — > 0. This assumption fails for the perturbed two di- 
mensional plane-wave Riemann problem. Consider a discontinuity line S. Without 
restriction of generality we may assume that for a sufficiently small portion, the 
discontinuity line S is perpendicular to the x-axis and we may assume that the 
flow is smooth in the y-direction. Furthermore for a sufficiently small time interval 
the shock speed s may be considered as constant. Let x = s(t,y) = s(y)t be the 
spatial-time discontinuity surface across which u has a jump. We obtain from the 
conservation law ([8]) for A x > 

g r+A x /2 ^ t)dx = _ [f ( u (A,/2, y, t)) - f (u(-A a /2, y, t))} 

Ot J-A x /2 

r+A x /2 Q 

J-A x /2 dy 
The first integral can be rewritten as 

Q ,-+A x /2 r+A x /2 Q 

— u(x,y,t) dx = -s[u(s+,y,t) - u(s-,y,t)] + / — u(x, y, t n ) dx 

Ot J - Ax /2 J-A x /2 Ot 

where s+ and s— denote one-sided limits at the discontinuity line. Thus if (|47j) 
vanishes for A x — > we obtain the jump conditions <\15h from the last two equa- 
tions. The integral does not vanish, if density perturbations can grow infinitely 
fast in time. 



Viscosity and heat conduction in the Navier-Stokes equations strongly affects steep 
gradients and perturbations. In reality a contact surface cannot be maintained for 
an appreciable length of time; (viscosity and) heat conduction between the per- 
manently adjacent particles on either side of the discontinuity would soon make 
the idealized assumption unrealistic. While gas particles crossing a shock front are 
exposed to heat conduction for only a very short time, those that remain adjacent 
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on either side of a contact surface are exposed to heat conduction all the time. 
Hence a contact layer will gradually fade out; see [CF1948]. Therefore some odd 
behavior of the solution must be expected for a mathematical idealized contact 
discontinuity. 



3 Viscosity and Heat Conduction 

It is well known that viscosity and heat conduction can attenuate small scale 
oscillations. In this section we study the effect of heat conduction and viscosity 
on the solutions (|26|) derived in the previous section. Including viscosity in our 
considerations we obtain from the Navier-Stokes equation for a divergence free 
flow 

d d d , A . 

dt P + U d- X P + V d-y P = ° ^ 
d d d Id fx d 2 d 2 

u + u in u + v ir u = -zinP + -Ain2 u +i^ u ) ( 48b ) 



dt dx dy pdx p dx 2 dy 2 

d d d Id p, , d 2 d 2 

—v + u—v + v—v = — — p + -{-^v + — - 
ot ox uy p ay p ox z ay 



d d d d 2 d 2 

c ^m T + u l>i T + v Ty T] = k ^ T + W T) (48d) 

< 48e » 

see e.g. [ATP1984;Section 5]. T is the temperature given by ©. For a constant 
x-component of the velocity u = uq we obtain from the divergence condition (|48cl) , 
that the tangential component v does not depend on y and equations ([15]) reduce 
to 

¥t p + U0 ^ p + W = ° (49a) 

d d 1 d p d 2 

—v + uq—v = 7rP+-7T^ v ( 49b ) 

ot ox p oy p Ox A 

|p + _ ( 7 - + ( 7 - m ^ 2 T + *,T) (49c) 
ox oy 
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where we used (jU) and ([5]) and assumed that the viscosity p, and the coefficient of 
thermal conductivity k are positive constants. We obtain for a density of the form 
(|26ap with u r = uq and v 6 = v for x < u$t with ()49bp 

d d d 

, d d d d 

at ax ay ay 

_ .1 d n d 2 . d 
p dy p dx 2 dy 

For t > (|49a|) requires 

|-f-W = or |°=0 (50) 
ay ax" 1 077 

Since the x-component of the velocity is constant, the pressure cannot depend on 
x. Since v only depends on the spatial coordinate x, this condition can only be 
satisfied if 

dy- p = ^ V = m OT (51) 

for a function f(t). Therefore either the density does not depend on y or v is a 
quadratic function in x; i.e. for a spatial oscillatory tangential velocity v, the den- 
sity p cannot vary in the y-direction. If viscosity is taken into account, perturbed 
solutions of the form ([25]) with an oscillatory tangential velocity Vq are excluded. 
In a constant pressure flow any density perturbation is equivalent to an entropy 
perturbation. Since is for x ^ u r t a smooth solutions of the Euler equations, 
these solutions cannot be excluded through the established entropy conditions for 
weak solutions. Thus we may regard (|51|) as an additional "plane wave condition" 
for the Euler equations. A density perturbation of the form ([26]) which violates 
(|5ip will be denoted as an "inviscid entropy perturbation". 
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For a constant pressure p = po equations (USD reduce to 



d d d 

d d d p . d 2 

—u + u—u + v—u 
at ox dy 

d d d 



u _j_ 9 2 ^ 

p dx 2 dy 2 



9 + ° + 9 _ p,d 2 + 9 2 
dt dx dy p dx 2 dy 2 

p d 



9 2 T + 9 2 ^ _ 2^/5 

d d 
—u + —v = 
ox oy 



•f _ 2^(—v) 2 - ~(—v + —u) 
k y dx^ k dy k dx dy 



(52a) 
(52b) 
(52c) 
(52d) 
(52e) 



where we used (jl]) to rewrite the temperature as a function of the pressure and 
the density and used (|52ap to obtain (|52dj) . For a given velocity field equation 
(|52d|) is a Poisson equation for the temperature; where the ratio of the viscosity 
and thermal conductivity in (I52dj) can be computed from the Prandtl number Pr 
and the specific heat at constant pressure c p by 

f = - < 53 > 

The ratio Pr/c p is approximately constant for most gases. For air at standard 
conditions Pr/c p ~ 0,00072. 

The temperature and constant pressure defines the density via the equation of 
state; i.e. 

For a non vanishing coefficient of thermal conductivity k ()52dj) imposes a regularity 
condition on the density, which is not present in an flow governed by the Euler 
equations ([8]); e.g. using the Weyl Lemma [WAl994;Section 9] we obtain from 
(|52d|) for a velocity field which is constant outside a bounded set, that the density 
is as smooth as the velocity field, whereas the Euler equations admits contact 
discontinuities for such a velocity field; see Proposition 1. 

A similiar regularity condition holds for the pressure in the incompressible Navier- 
Stokes equations. For a constant density p = po the equations (|48|) reduce to 
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(55d) 



where we used (|54p and (JSJ to express the temperature as a function of the den- 
sity and pressure. Equation (|55a|) . (|55bp and (|55dp constitute the incompressible 
Navier-Stokes equations. The pressure in these equations can be computed from 
the velocity field through a Poisson equation 



see e.g. [KL1989; Section 9.1.3]. Again using the Weyl Lemma we obtain from 
(|56|) for a velocity field which is constant outside a bounded set, that the pressure 
is as smooth as the velocity field. Inserting (|56p into ()55c[) we obtain for a given 
smooth velocity field a inhomogeneous advection equation for the pressure. This 
equation defines the time evolution of the pressure and is generally neglected in 
an incompressible flow. Consider a constant pressure, constant density flow with a 
constant x-component of the velocity u = u r . We obtain from (|56|) for a constant 
pressure, that the functional determinate of (u,v) vanishes. Therefore v = <&(«) 
for some smooth function and we obtain that the tangential velocity is constant 
as well. From proposition 1 follows that (u r ,v) with v = &(x — u r t) is a velocity 
field for the Euler equations with a constant pressure and constant density. These 
solutions are excluded for a finite viscosity fi and coefficient of thermal conductiv- 
ity k in the divergence free Navier-Stokes equations (|48l ) . 

Numerical methods for the Euler equations employ an artificial numerical viscosity 
model to resolve discontinuities. If this numerical viscosity model is not properly 
related to the physical viscosity and thermal conductivity, numerical artefact 's 
may be introduced into a approximate solution, which are not related to a real 
viscous, heat conducting flow. 




(56) 
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4 Instability of Godunov's method 



In a first order finite volume methods the assumption is generally made, that the 
solution is constant inside a cell at a time-level t = t n ; e.g. 

u(x,i) := Uij for x e [zi_i/ 2 , Xi+1/2] x [yi-1/2, 2/1+1/2] (57) 

and it is assumed that discontinuities are moved to the cell boundary. We assume 
for simplicity that Aj = [^-1/2)^1+1/2] x [Hi- 1/2 > 2/1+1/2] is a rectangle defined 
through a constant cartesian grid x% = iA x and yj = jA y and denote by 

the cell-average, where K(-Djj) is the Volume of -Djj. Finite- Volume Godunov- 
type methods are derived from the integral form of the conservation law. The 
integral form (|13j) can be rewritten as 

4** = -TTTtVt / H (0 ■ F(u R (x(0,t+; H@)) d£ (59) 

at V \ D i,j) JdD it j 

This equations says that we can evolve the cell averages in time, by solving one- 
dimensional Riemann problems at a the cell boundary at time t = t n and then 
solve the system of ordinary differential equation (|59[) to obtain the cell average 
at time t = t n + r, r > 0. Taking for granted that the solution of the Riemann 
problem at the cell interface can be locally advanced in time; i.e. we require that 

^-u R (m,t;n(0)\m (60) 
should exist at the cell-boundary dDij. 

In Godunov's method a plane- wave Riemann problem is solved at the cell bound- 
ary, for example at the cell-boundary (xi + i/2, Uj) at time t = t n in the x-direction. 
Denote by u R (x i+ i/ 2 , Vj, t) the solution of the plane-wave Riemann problem in 
the x-direction at the cell boundary, then u™ +1/ , 2j . := lim t ^t n + u" R (x i+1 / 2 , Vj, t) is 
computed and used to evaluate the physical flux function. 

In the last two sections we saw that a constant pressure region which a contact dis- 
continuity is unstable under perturbations. In this section we analyze Godunov's 
method for this critical region. Assume that the velocity u and the pressure p is 
constant, then the solution of the Riemann problem in the x-direction consists of 
contact discontinuity and we obtain 

p/ „s I Uj+i ,• for u < . . 

u R (x i+1/2 , yj ,t n ) = \ (61) 
u,- 1 for (J < u 
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If v depends only on the spatial coordinate x and p is constant, then the solution 
in the y-direction is also a contact discontinuity and we obtain 

u {Xi,y j+1/2 ,t ) - < _ (62) 



Uij for < Vi 



Therefore for u < and Vj < we obtain 



which is equivalent to the four difference equations 

7" T" 



Since the second equation is just the first equation multiplied with the constant 
velocity u, this can be reduced to three difference equations 

pit 1 = Pi,j - i~ u iPUu - Kj) - i-i>i<tfj+i - pIj) 



Eft 1 ~ E tj ~ -£~ u ( E i+i,j ~ E i,j) ~ -^Vi( E Zj+i ~ E i,j 

x y 



This is a discrete approximation to 



Ft P + U ^ P + V Yy P = ° (64a) 
d d 

£* + »s* +e £* =0 (64c) 



Rewriting the second equation as 



p[ ¥t v + u Tx v] + v[ ¥t p + u T x p + W ] = 
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we obtain that (|64p . is equivalent to (|24|) in this flow region. This holds also for 
u < 0, v > 0, u > 0, v < and u > 0, v > 0. Therefore Godunov's scheme is a 
discrete approximation to (|24|) in a regions with a constant pressure p, constant 
x-component of the velocity u and a y-component of the velocity v, which does 
only depend on the spatial coordinate x. 

For a stationary shock wave at x-^i^ = with a shock line perpendicular to the 
x-axis, the x-component of the flux function of Godunov's scheme satisfies 

ff /2)i = f(u^.) = f(ug^) (65) 

i.e. one-dimensional stationary shocks are resolved exactly. Therefore Godunov's 
method is a discrete approximation to (|24l) behind the shock with boundary con- 
ditions 

p(0, t) = const. 
v(Q, t) = const. 

where u r in (|24|) is the constant normal component of the velocity behind the 
shock. If the shock is near stationary these boundary conditions introduce per- 
turbation. Thus on a sufficiently fine grid entropy perturbations of the form (|26|) 
and acoustic perturbations of the form (I45p are resolved by Godunov's method 
behind a near stationary shock. In the last section we saw that inviscid entropy 
perturbations do not relate to a real viscous flow. Flow structures not related to 
a real viscous, heat conducting flow may therefore appear in numerical solutions 
of Godunov's method. 

Numerical examples for a plane shock wave aligned with the grid, which is moving 
down a duct are given in [QUI 1994; Figure 5]. At the grid center line, a small per- 
turbation is introduced in the computation. Downstream of the shock an unstable 
density profile develops, which over time leads to an unstable numerical shock 
front. If we associate the center of the duct with the x-axis, then the perturbation 
depend behind the shock front on y. This numerical example reflects the situation 
discussed analytically at the beginning of section 2. 

Roe's method [Roel980] is more prone to generate these perturbations then Go- 
dunov's method, due to the fact that a rarefaction wave is replaced by a rarefaction 
shock. Therefore noise from the rarefaction shock, can also lead to an unstable 
growth of the density, in addition to the noise from the shock. 
An example from an aerodynamic simulation, which results in incorrect numerical 
results is given in [PEIM1988], for a bow shock over a blunt body placed in a high 
Mach number flow. Along the stagnation line the bow shock is approximately 
aligned with grid used for the calculation. A perturbation of the shock profile is 
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given through the curvature of the shock. At the stagnation point we have approx- 
imately a plane-wave near stationary shock wave, with a disturbed shock profile. 
This again is the situation discussed at the beginning of section 2. 
Based on the numerical observation, that shock capturing methods which try to 
capture contact discontinuities exactly, generally suffer from failings, a link be- 
tween the carbuncle phenomenon and the resolution of the contact discontinuities 
was suggested by Gressier and Moschetta [GRE1998]. 

The dissipation model for Godunov's scheme was studied by Xu [XU1999] in a 
series of numerical experiments. He concluded that Godunov's method gives accu- 
rate results in both unsteady shock structure and boundary layer calculations, but 
that the absence of dissipation in the gas evolution model in Godunov's scheme 
amplifies post-shock oscillations. He found that the numerical dissipation model 
for Godunov's method is in the multidimensional case mesh-oriented and not con- 
sistent with the Navier-Stokes equations. Xu also mentioned that it is well known 
that the inviscid Euler equations cannot give a correct representation of the fluid 
motion in the discontinuity flow region. Which is consistent with our analysis re- 
garding a contact discontinuity line (vortex sheet). 

Real weak shock fronts are transition layers of finite width and the representation 
of weak shock fronts through a discontinuity line in the Euler equations is a mathe- 
matical approximation. A justification for this approximation was given by Majda 
[MA1983; Proposition 3]. Based on a linear stability analysis Majda found that 
planar compressive shock fronts in an ideal gas are uniformly stable. In [LL1959; 
§87] it is furthermore noted that the front depth of a non weak shocks is so small, 
that a transition layer becomes meaningless. Also by no means obvious, it may 
be safely assumed, that the mathematical representation of a single smooth shock 
front in an ideal gas through a discontinuity line resp. surface is reasonable. Since 
Majda 's stability analysis does not apply to a contact discontinuity line/surface 
and we assumed in this paper that shock fronts only generate perturbations with- 
out assuming an unstable growth at the smooth shock front itself, this result is 
consistent with the analysis in this paper. 

The source of small perturbations in Godunovs methods is the displacement of the 
shock curves at the cell-boundary and the nonlinear interaction of the dependent 
variables in the numerical shock layer. A non stationary or stationary displaced 
shock wave is approximated through a smeared profile with at least one intermedi- 
ate cell. Whenever the smeared shock profile changes, perturbations are generated 
from the characteristic fields. If the shock curve is not exactly a plane- wave in the 
x-direction, the shock-curvature will introduce an y dependence in the perturba- 
tions. In regions of low (numerical) viscosity, the instability under perturbation of 
the constant pressure region behind the shock will result in an unstable flow - if a 
contact discontinuity is present. 
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For the case reported by Perry and Imlay [PEIM1988] entropy perturbations are 
generated a the shock. Since according to the proposition in section 2, entropy 
perturbations can only exits if the tangential velocity does not depend on y, these 
perturbations are initially restricted to the symmetry axis. If the flow is near 
stationary, cell-averaging of the form (|58p in the projection state in Godunov's 
method leads to density perturbation in these cells. At the boundary of these 
cells Godunov's scheme then resolves a contact discontinuity, which results in an 
unstable flow according to the analysis in section 2 and 4. The unstable flow be- 
hind the shock interacts with the near stationary shock front resulting in carbuncle 
structures. Since the interaction destroy 's also the smoothness of the shock front 
itself, the stability analysis of Majda no longer applies. 

For the case reported by Quirk [QUI 1994; Figure 5] perturbation are introduced 
externally at the grid center line. The grid center line corresponds to the symme- 
try line in results of [PEIM1988]. In [QUI 1994; Figure 5a] perturbations are first 
visible behind the plane shock, but no carbuncle is visible at this stage. Only after 
the shock has propagated further down the duct the carbuncle becomes visible in 
front of the shock. 

If discontinuous density or velocity perturbation are artificially introduced in a 
background flow then according to the analysis in section 2, unstable structures 
can even be generated in numerical solutions of the isentropic Euler equations. 
This was demonstrated by Elling [VE2006] for Godunov's method. In his case a 
steady plane shock line parallel to the y-axis is perturbed through a one-cell-high 
filament along the x-axis in front of the shock. In the filament the normal com- 
ponent of the velocity u = u is set to zero. The flow in the filament is therefore a 
constant pressure flow governed by (|32h with a constant normal velocity compo- 
nent v = v. From proposition 1 follows (with the tangential and normal velocity 
interchanged) that u = u(y,t). The linear stability analysis results again in an 
advection equation with a singular source term of the form 

d _ d 

—u + v—u = -v(u r -u l )5{y-vt) (66) 
at ay 

This shows that the flow in the filament is unstable. Interactions of this unsta- 
ble flow with the shock wave results in a carbuncle instability. Elling noted that 
the carbuncle instability can also be observed for the local Lax-Freidrich/Rusanov 
and Osher-Solomon schemes. His conjecture is that carbuncles can be related to a 
special class of non-physical entropy solutions for the continuum equations, which 
is supported by our analysis. 

A discussion of the carbuncle instability for several numerical methods, can be 
found in [DMG;2004] and [Roe2007]. The focus in these papers is on the pertur- 
bations generated through a numerical shock profile. The cause for the carbuncle 
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instability is often related to the numerical shock profile itself and the instability 
also denoted a as a shock instability. However, in this paper it is shown, that the 
instability may be related to the central constant pressure region in the plane wave 
Riemann problem. This is an intrinsic instability of the Euler equations and not 
related to a particular numerical method. The numerical shock profile manifests 
this instability in regions of low numerical dissipation through the generation of 
perturbations. This manifestation of the instability, is problem and method de- 
pendent. 



5 Summary and Conclusion 

In this paper it is proved that a constant pressure flow region, governed by the 
hyperbolic conservations laws for an ideal gas, is unstable under perturbations if a 
discontinuity is present. The instability is immanent to the linear degenerate fields 
in the multidimensional hyperbolic conservation laws, with a double eigenvalue. 
The mathematical idealized assumption of a zero thickness contact discontinuity 
in the plane wave Riemann problem is unrealistic for real flows. If infinitesimal 
viscosity and thermal conductivity are taken into account, inviscid entropy per- 
turbations are excluded for an oscillatory tangential velocity in a perturbed plane 
wave Riemann problem. Since smooth entropy disturbances in a constant pres- 
sure flow are transported with the fluid, they cannot be excluded through the 
established entropy conditions for weak solutions of hyperbolic conservations laws. 
Additional conditions are required to guarantee that a solution of the plane wave 
Riemann problem is as a limit solution of the Navier-Stokes equations for a van- 
ishing viscosity and thermal conductivity. 

Immanent perturbations (acoustic and entropy) generated at a perturbed shock 
line and the numerical impossibility to differentiate exactly between admissible 
and non admissible entropy perturbations in a high Reynols number (low vis- 
cous) flow, are challenges for discontinuity resolving methods. Godunov's method 
closely relates to the physics of the Euler equations. Numerical artefact 's observed 
in numerical solutions for Godunov's method are a numerical manifestation of the 
immanent instability in the Euler equations. The HLLE scheme [EIN1988] on 
the other hand can be regarded as a vortex sheet averaged approximation to the 
Navier-Stokes equations for high Reynolds number flows. 
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